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Abstract 

We use a two-state ratchet model to study the cooperative bidirectional motion of 
molecular motors on cytoskeletal tracks with randomly alternating polarities. Our model is 
based on a previously proposed model [Badoual et al., Proc. Natl. Acad. Sci. USA 99, 6696 
(2002)] for collective motor dynamics and, in addition, takes into account the cooperativity 
effect arising from the elastic tension that develops in the cytoskeletal track due to the joint 
action of the walking motors. We show, both computationally and analytically, that this 
additional cooperativity effect leads to a dramatic reduction in the characteristic reversal 
time of the bidirectional motion, especially in systems with a large number of motors. We 
also find that bidirectional motion takes place only on (almost) a-polar tracks, while on 
even slightly polar tracks the motion is unidirectional. We argue that the origin of these 
observations is the sensitive dependence of the cooperative dynamics on the difference 
between the number of motors typically working in and against the instantaneous direction 
of motion. 
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1 Introduction 



Many cellular processes such as cell motility and mitosis require the cooperative work of many 
motors in order to preserve continuous motion and force generation [1]. Muscle contraction, for 
example, involves the simultaneous action of hundreds of myosin 11 motors pulling on attached 
actin filaments and causing them to slide against each other [2]. Groups of myosin 11 motors 
are also responsible for the contraction of the contractile ring during cytokinesis [3]. In certain 
biological systems, cooperative behavior of molecular motors produces oscillatory motion. In 
some insects, for instance, cooperative behavior of molecular motors leads to oscillations of the 
fiight muscles [4]. Another example is the oscillatory motion of axonemal cilia and flagella, 
which is believed to be generated by a large number of interacting dynein motors [5,6]. Finally, 
cooperative action of motors is required for the extraction of membrane tubes from vesicles [7, 8] 
One of the more interesting outcomes of cooperative action of molecular motors is their ability 
to induce bidirectional motion. "Back and forth" dynamic has been observed in various motility 
assays including: (i) myosin 11 motors walking on actin tracks with randomly alternating po- 
larities [9], (ii) NKll (kinesin related Ned mutants which individually exhibit random motion 
with no preferred directionality) moving on microtubules (MTs) [10], (iii) mixed population of 
plus-end (kinesin-5 KLP61F) and minus-end (Ned) driven motors acting on MTs [11], and (iv) 
myosin 11 motors walking on actin filaments in the presence of external stalling forces [12]. Re- 
versible transport of organelles through the combined action of kinesin 11, dynein, and myosin 
V has been also observed in Xenopus melanophores [13]. In the latter example, the kinesin and 
dynein move the organelle in opposite directions along MTs, while the myosin motors (which 
take the organelle on occasional "detours" along the actin filaments) function as "molecular 
ratchets" , controlling the directionality of the movement along the MT transport system. From 
a theoretical point of view, cooperative dynamics of molecular motors and, in particular, bidi- 
rectional movement, have been investigated using several distinct models. These models include: 
(i) lattice and continuum asymmetric exclusion models [14-20] , (ii) ratchet models of interacting 
particles moving in the presence of a periodic potential [16,21-24], and (iii) the tug-of-war model 
which has been recently proposed for describing the transport of cargo by the action of a few 
motors [25-28]. The common theme in these experimental and theoretical studies is the asso- 
ciation of bidirectionality with the competition between two populations of motors that work 
against each other to drive the system in opposite directions. The occasional reversals of the 
transport direction reflect the "victory" of one group over the other during the respective time 
intervals. The balance of power is shifting between the two motor parties as a result of stochas- 
tic events of binding and unbinding of motors to the cytoskeletal track. Without going into the 
details of the various existing models of cooperative bidirectional motion, we note that most of 
them assume that the motors interact mechanically but act independently, i.e., their binding to 
and unbinding from the track are uncorrelated. By further assuming that the attachment and 
detachment events of individual motors are Markovian, the distribution of "reversal times" (i.e., 
the durations of unidirectional intervals of motion) can be shown to take an exponential form 

p{6t) ^ eX.p{-St/Tre^r): (l) 

where Trev is the characteristic reversal time of the bidirectional motion. 

The magnitude of Trev can be taken as a measure for the degree of cooperativity between 
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the motors. The more cooperative the motors are, the more persistent is the movement and the 
longer are the periods of unidirectional transport. The run lengths (in each direction) of highly 
cooperative motors may be of a few microns even for non-processive motors like myosin II [9]. 
As noted above, the mechanical coupling between the motors is sufficient for the generation of 
highly cooperative bidirectional motion, even if the motors attach to/detach from the track in 
an uncorrelated fashion. This has been demonstrated theoretically by Badoual et al. some years 
ago [23]. A slightly modified version of this model is illustrated schematically in Fig. [H The 
model considers the one-dimensional motion of a group of point particles (representing the 
motors) connected (mechanically coupled) to a rigid rod with equal spacing q. The cytoskeletal 
track is represented by a periodic saw-tooth potential, U{x), with period / and height H. The 
model requires that q is larger than and incommensurate with I. The motors are identical and 
walk on a track which is globally a-polar and, thus, does not permit net transport to the right 
or left over large time scales. The temporal direction of motion is determined by the net force 
generated by all the motors. The local polarity of the track is represented by an additional force 
of size /ran (dcuotcd by a horizontal arrow in each periodic unit in Fig. [T]) which, within each 
unit of the periodic potential, points to the right or to the left. The globally a-polar nature of 
the track is ensured by requiring that the sum of these random forces vanishes. 

The instantaneous force between the track and the motors is given by the sum of all the 
forces acting on the individual motors: 



1=1 1=1 



dx 



Q{t), (2) 



where Xi = Xi + [i — l)q is the coordinate of the i-th motor. The two terms in the square brackets 
represent the forces due to the symmetric saw-tooth potential and the additional random local 
forces acting in each periodic unit. The function Ci(t) takes two possible values, or 1, depending 
on whether the motor i is detached from or attached to the track, respectively, at time t. The 
motors change their binding states (0 - detached; 1 - attached) independently of each other, 
according to the following rules: We define an interval of size 2a < I centered around the 
potential minima (the gray shaded area in Fig. [1]). If located in one of these regions, an attached 
motor may become detached (1 — 0) with a probability per unit time Ui. Conversely, a detached 
motor may attach to the track (0 — > 1) with transition rate uj2 only if located outside this region 
of size 2a. 

At each instance, the group velocity of the motors is proportional to the total force exerted 
by the motors (Eq. [2]) 

vit) = Ftot(t)/A, (3) 

where the friction coefficient. A, depends mainly on motors attached to the track and is, there- 
fore, taken proportional to the number of connected motors, A'c < A^ at time t: \ = XqNc- 
Because the track is globally a-polar, it is clear that the motors exhibit "back and forth" motion 
with vanishing mean velocity and displacement. The characteristic time of movement in each 
direction, r^ev, may nevertheless be macroscopically large. The origin of this feature (which re- 
fiects the cooperative character of the motors' action) can be explained as follows: The stochastic 
equations of motion of our model system have two identical (except for sign reversal) steady state 
solutions corresponding to right and left motion of the mechanically coupled motors. Each of 
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these solutions is characterized by A^^c = ■ P (P < 1) connected motors. The connected motors 
are partitioned to and = Nc — < motors that, respectively, support and object 
the motion. Let us define the excess number of motors working in the direction of the motion as 
N -A = iV+ — iV_, where A will be termed the "bias parameter" . Notice that P and A denote the 
averages of quantities (which we, respectively, denote by P(t) and A(t)) whose values fiuctuate 
in time due to the stochastic binding and unbinding of motors. To switch the direction of mo- 
tion, A{t) must vanish, and the occurrence probability of this event, n(A(t) = 0) ~ (Trev)~^- In 
section [3] we derive an approximate expression for n(A(t) = 0) and show that the mean reversal 
time of the bidirectional motion increases exponentially with the size of the system: 



Thus, for sufficiently large N, the characteristic reversal time of the bidirectional motion becomes 
macroscopically large. In the "thermodynamic limit" (A^ oo), Trev diverges and the motion 
persists in the direction chosen at random at the initial time. 

The validity of Eq. H] was recently tested using a motility assay in which myosin II motors 
drive the motion of globally a-polar actin bundles [9]. In contrast to the predicted exponential 
dependence of Xrev on the number of working motors, the experimentally measured reversal times 
showed no dependence on A^. The apparent disagreement between the theoretical model and the 
experimental results can be reconciled by noting that Eq. H] describes exponential growth of r^ev 
with A^ only when P and A are themselves independent of A^. This is indeed the case in the 
original model presented by Badoual et al. [23], where both the on {1JJ2) and off (tui) rates do not 
depend on A^. In ref. [9] we introduced a slightly modified version of Badoual's model, which to 
a large extent explains the experimentally observed independence of Trev on A^. We argued that 
the origin of this behavior can be attributed to the tension developed in the actin track due to 
the action of the attached myosin II motors. An increase in the number of attached motors leads 
to an increase in the mechanical load which, in turn, leads to an increase in the detachment 
rate of the motors, as already suggested in models of muscle contraction [29-32]. But unlike 
most previous studies where the myosin conformational energy was calculated, in ref. [9] we 
considered the elastic energy stored in the actin track. Within a mean field approximation, this 
energy scales as ~ {F^^^)/ksp, where ksp is the effective spring constant of the track and Ftot is 
the total force exerted by the motors (see Eq. [3]). The total force is the sum of A'c random forces 
working in opposite directions. Therefore, the mean force (-Ftot) = 0, while (-Ftot) scales linearly 
with Nc- The spring constant is inversely proportional to the length of the track, i.e. to the size 
of the system and to the total number of motors A^. We thus conclude that the mean elastic 
energy of the actin scales like i^/fc^T ~ NNc, which means that the detachment of a motor 
[Nc Nc — 1) leads, on average, to an energy gain dE/ksT = —aN [a is some dimensionless 
number). This effect can be incorporated within the model by introducing an additional off rate, 

= tUg exp(aA^), outside the gray shaded area in Fig. 1 (i.e., around the potential maxima). 
Thus, in the modified model, the on and off rates within each unit cell {—1/2 < x < +//2) of 
the periodic potential are given by 



rev 



1-P + VP^ - A2 



(4) 




(5) 
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and 

f . j uji \x\ <a , . 

y u;3exp(a;A') a < |s| < 1/2. 

The dependence of uJos{x) on is another, indirect, manifestation of cooperativity between the 
motors which is mediated through the forces that the motors jointly exert on the actin track. 
Notice, however, that the cooperative action of the motors does not lead (within the modified 
model) to correlations between attachment and detachment events of different motors. Moreover, 
the attachment/detachment rates do not depend on the number of connected motors and, thus, 
are fixed over time. In ref. [9] we compared the predictions of the original and modified models 
for model parameters corresponding to the myosin Il-actin motility assay. In the former, the 
reversal times grew exponentially with N from Trev ~ 1 sec for N ~ 1000, to Trev ~ 10'^ sec 
for ~ 3000. The modified model showed a much better agreement with the experimental 
results [9]. Specifically, the reversal time did not grow exponentially with N but rather showed 
a weak maximum around ~ 2000, where Trev ~ 10 sec. 

In this paper we present a more detailed account of the theoretical model. Several aspects 
of the model not studied in ref. [9] will be discussed in section [2], including the dynamics on 
non-random and slightly polar tracks. The steady state solutions of the bidirectional motion 
are derived analytically in section [31 We use these solutions to estimate the reversal times and 
compare our analytical predictions with the computational results. We summarize and discuss 
the results in section HI 



2 Computational results 

The model has been presented in details in ref. [9] , as well as above in section [H The model 
parameters used in this paper are summarized in Table [H The choice of these values which 
represent various chemical and physical parameters of the myosin Il-actin system is explained 
in detail in ref. [9] . The model features two new parameters which do not appear in the original 
model [23]. The off rate represents the probability of a motor to detach from the track 
without completing a unit step, and its value was estimated in ref. [9] by noting that in the 
absence of an elastic load, the probability of such an event is 1-2 orders of magnitude smaller 
than the complementary probability that the attached motor will execute the step. The constant 
a depends on the effective elastic spring constant of the basic actin unit (monomer) as well as 
on the magnitude of the forces that the motors typically apply on the track (see Eq. 5 in 
ref. [9]). The best fit to the experimental data is achieved with a = 0.0018 [9], which gives a 
weak non-monotonic dependence of Trev on A^. Here we set a = 0.002, which gives somewhat 
poorer agreement with the experimental results but which highlights the differences between 
Badoual's model [23] and the newly proposed model that takes the elasticity of the actin track 
into account. As noted above, the consideration of the elastic properties of the actin and the 
cooperative nature of the action of the working motors considerably improves the results of the 
original model by replacing the very strong exponential dependence of Trev on A^ with a much 
weaker non-monotonic dependence. 

To simulate conditions corresponding to dynamics on a-polar tracks, we randomly chose the 
direction of random force (the force representing the local polarity of the track, see horizontal 
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parameter 



value 



/ - period length of the potential 

q - spacing between adjacent motors 

2a - size of the gray shaded area (see Fig. 7A) 

2H/1 - force due to the periodic potential 

/ran " raudom force in each unit cell 

(cui)-i (see Eq. El 

(cua)-^ (see Eq. E) 

(cu^)-^ (see Eq. ED 

a (see Eq. [6]) 

Aq - friction coefficient per connected motor 



0.5 msec 
33 msec 
7500 msec 
0.002 

85 X 10^ kg/sec 



(57r/12)/ ~ 6.54 nm 
3.8 nm 



5 pN 
1 pN 



5 nm 



Table 1: Values of the model parameters as used in our simulations. 



arrows in Fig. [T]) in each unit cell, but discarded the tracks at which the sum of random forces 
did not exactly vanish. We computationally measured characteristic reversal time, 
function of in the range of 400 < N < 2400. For each value of A^, we generated 40 different 
realizations of random tracks (each of which consisting of M ~ {q/l)N units with periodic 
boundary conditions) and simulated the associated dynamics for a total period of 2 ■ 10^ seconds. 
During this period of time we followed the changes in the direction of motion and calculated the 
probability distribution function (PDF) of the reversal times. The characteristic reversal time 
corresponding to each random track was extracted by fitting the PDF to an exponential form 
(see Eq. [1]), as demonstrated in Fig. [2]A. Fig. [2]B summarizes our results, where here for each N 
the reversal time plotted (denoted by (r^ev)) is the average of Tj-^v calculated for the different track 
realizations. The error bars represent the standard deviation of Trev between realizations. The 
data points depicted in solid circles correspond to a = 0.002, while the open circles correspond 
to a = 0, i.e., to the model originally presented in ref. [23] where the on and off rates defined in 
Eqs. [5] and [6] are independent of N. As predicted by Eq. H] and indicated by the straight line in 
Fig. 2B, for a = the mean reversal time (rrev) exhibits a very strong exponential dependence 
on A^. Because of this very rapid increase of (TVev) with A^, the reversal times (in the a = 
case) could not be accurately measured for A^ > 1800 . Based on the exponential fit (solid line 
in Fig. 2B), we estimate that for A^ = 2400 the mean reversal time will be of the order of a 
few hours. In contrast, the calculated (ttcv) corresponding to a = 0.002 show a non-monotonic 
dependence on A^. The computed (rrev) are much smaller in this case, and fall below 1 minute 
for all values of A^. In ref. [9], we used slightly different values of the model parameters than 
those given in Tabled! For the model parameters in ref. [9], the variation of (rrev) with A^ was 
even weaker than in Fig. 2B and the largest computed (rrev) < 12 sec. These computational 
results were in a very good quantitative agreement with the experimental results of the in vivo 
actin-myosin motility assay. 

Our simulations reveal surprisingly large variations in the rrev values between random tracks 
of similar size (see error bars in Fig. [2]B). The origin of these variations lies in the fact that the 
spacing between motors is larger than the periodicity of the ratchet potential {q > I) and, thus, 
only A^ out of M ^ (g//)A^ unit cells are "occupied" with motors (which may be either connected 
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or disconnected) at each instance. Thus, although the track is perfectly a-polar and contains an 
equal number of cells with random forces pointing in both directions, the subset of occupied cells 
may have net polarity which constantly changes with time as the motors move collectively along 
the track. The direction of the net polarity of the occupied cells is also the instantaneous preferred 
direction of motion. Therefore, the temporal variations in the net polarity must be correlated 
with the changes in the directionality of the motion. We thus expect tracks on which the net 
polarity changes more frequently to have smaller Tj-ev The effect of net polarity fluctuations 
does not occur when the motion takes place on periodic a-polar tracks, because in this case the 
equally spaced motors occupy equal number of cells with left- and right-pointing random forces. 
Therefore, the reversal times on periodic a-polar tracks are expected to be (i) independent of 
the periodicity of the track, and (ii) larger than the reversal time on random a-polar tracks. 
These predictions are fully corroborated by the results from simulations with two very distinct 
periodic tracks - one with period 2 (i.e., where the local random force changes its sign every 
unit cell) and one with period M (i.e., when the track is divided into two equal domains of 
opposite polarities). The results from these two sets of simulations are denoted by triangles in 
Fig. [2]B. The reversal times of both periodic tracks are nearly indistinguishable from each other 
(the differences between them are smaller than the size of the symbols) and are larger than the 
reversal times measured for all the random tracks of similar size. 

What happens when the simulated track is not perfectly a-polar and the number unit cells 
in which the random force is pointing in one direction is slightly larger than in the opposite 
direction? Obviously, the nature of motion is expected to gradually change from bidirectional 
to unidirectional. In order to investigate this transition between two types of dynamics, we 
simulated the dynamics of motors on tracks in which the fraction of random forces pointing in 
one direction, pi, is slightly larger than 0.5. In the simulations, we fixed the number of motors to 
= 1000 and varied the difference D = pi — {1 — pi) = 2pi — 1 between the fractions of random 
forces pointing in the favored and disfavored directions. For each track, the simulation data 
was analyzed in a manner similar to that described above for a-polar tracks, i.e., by fitting the 
PDF of time intervals of unidirectional motion to an exponential function. There is one notable 
difference, however, between the analysis of the results for a-polar and for slightly polar tracks. 
In the latter case, two PDFs, one corresponding for each direction of motion, must be generated 
with different characteristic reversal times. The motion in the preferred direction is characterized 
by the larger reversal time, Trev-i, which increases with D. Conversely, the smaller reversal time, 
Trev-s, Corresponding to the motion in the opposite unpreferred direction decreases with D. These 
observations are summarized in Figs. [3] and |H In Fig. [3l the PDFs corresponding to tracks with 
D values varying from to 0.05 are shown. In the a-polar case D = 0, the two PDFs coincide 
with each others, and the velocity histogram (see inset) is bimodal. As D increases, the two 
PDFs become increasingly different - the one corresponding to the preferred direction of motion 
becomes flatter (due to the increase in Trev-i), while the other one gets more peaked at the origin 
(as a result of the decrease in Tj-cv-s )• The fact that the motors spend larger time intervals 
moving in one direction is also reflected in the corresponding velocity histograms which become 
less and less symmetric. The results presented in FigJS] are obtained from simulations of six 
different track realizations, one for each different value of D. The mean reversal times, (rj-ev-i) 
and (Trev-s), obtained by averaging the reversal times computed for 8 track realizations for each 
value of D, are shown in Fig. HI For D = 0.05, (rrcv-i) ~ 10(ri.cv-s)- For even larger values of D, 
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the dynamics are essentially unidirectional, as intervals of motion in the unpreferred direction 
become very rare and short. We, thus, conclude that bidirectional motion can be observed only 
on a-polar or slightly polar tracks. This conclusion is directly related to the cooperativity of the 
motors which causes persistent motion that cannot be easily reversed. 



3 Analytical treatment 

In the following section we use mean field master equations to analyze the bidirectional motion 
exhibited by our computational model. The mean field approach corresponds to the limit N ^ 1 
where one can introduce the probability densities Pa,tt{x) and Pdct{x) of finding a motor in the 
attached or detached state, respectively, at position —1/2 < x < 1/2 within the unit cell of the 
periodic potential. These probability densities are the steady-state solutions of the following set 
of coupled master equations which govern the transitions between the two connectivity states: 

dtPs±t{x, t) + vd^ps±t{x, t) = -uJoa{x)ps.tt{x, t) + LUon{x)pdct{x, t) , . 

dtPdetix, t) + vd^PAct{x, t) = -UJon{x)pdct{x, t) + UJoa{x)p^ttix , t). 

In Eq. [71 uJon{,x) and ci;ofr(a;) denote the space-dependent on and off rates, and v is the group 
velocity of the motors. Because the spacing between the motors is incommensurate with the 
periodicity of the potential, the total spatial distribution is uniform in x for ^ 1: 

Patt {X,t)+ Pdct (X, t) = y . (8) 

Using Eq. [HI together with Eqs. [S] and [H] to define the on and off rates in Eq. [3 the following 

steady-state equation {dtP = 0) can be derived for Pa.u{x): 

j^^ dp^jx) ^ f -luip^ttix) for |a;| < a 

dx \ UJ2 - I {^2 + exp (aN)) p^ttix) for a < |a;| < //2. 

Eq [9] should be solved subject to the boundary condition that Patt(— V^) = Pa.tti}/'^) and the 
requirement that Pii,tt{x) is continuous anywhere in the interval —1/2 < x < 1/2, including at 
X = ±a. Several solutions are plotted in Fig. \5\ for 2a = 0.76/ (see Table [T]), = 0, and 
(u;i,ti;2) = {v/l,v/l) (thin sohd line), {5v/l,5v/l) (dashed line), and (30^//, 30f //) (thick sohd 
line). The solutions correspond to the case when the motors move to the right {v > 0) and, 
therefore, it is easy to understand why p^tt reaches its maximum at x = —a [just before the 
motors enter, from the left, into the central gray-shaded detachment interval {—a < x < a)] 
and its minimum at x = a [just before leaving the central detachment interval through the right 
side]. We also notice that when the off rate ui ^ v/l, patt drops very rapidly (exponentially) 
to near zero in the detachment interval. When the attachment rates 002 3> v/l, p^tt increases 
exponentially fast for x > a and rapidly reaches the maximum possible value patt = !//• The 
second steady state solution corresponding to the case when the motors move to the left {v < 0) 
is simply a mirror reflection of the first solution with respect to x = 0. 

The mean fraction of connected motors, P, can be obtained by integrating the function Patt (x) 
over the interval —1/2 < x < 1/2 

j.l/2 

P= / Patt(x)dx. (10) 

J -1/2 



8 



The population of connected motors can be divided into two groups: The connected motors 
which are located left to the minimum of the periodic potential (— //2 < x < 0) experience 
forces pushing them to the right, i.e., forces directed in their direction of motion. Conversely, 
attached motors which are located right to the minimum experience forces directed opposite to 
their direction of motion. Thus, the bias parameter A, previously defined as the excess mean 
fraction of motors supporting the motion, can be related to Patt by 



fO 1-1/2 
A/2 Jo 



'111 



In order to derive an expression for the reversal time of the dynamics, we now consider the 
fluctuations of the instantaneous bias parameter, A(t), around the mean value A. The motors 
may switch their direction of motion when A(t) = 0, i.e., when the motion momentarily stops. 
The occurrence probability of such an event can be related to the mean reversal by: 



Trev- [n(A (t) =0)]- 



(12) 



To estimate 11 (A (t) = 0) we proceed by noting that the probability of finding a motor attached 
left to the minimum of the potential, i.e. a motor experiencing a force directed in the direction of 
motion, is = (P+ A)/2. The probability that a motor is experiencing a force directed opposite 
to the direction of motion is P~ = (P — A)/2. The probability of having and < 
motors which, respectively, support and object to the motion can thus be approximated by the 
trinomial distribution function 



n{N+,N_) 



N+\NJ{N - N+ - N_)\ 



P + A 



P-A 



{N-N+~N^) 



The instantaneous bias is given by A{t) = {N^ — N_)/N, and the probability that A(t) 
be expressed as sum over the relevant terms in Eq. [13] for which = A^_|_ 



n(A(t) = 0) 



N/2 



N/2 

E 



'i\y{N-2t)\ 



P^- A- 
4 



;i-p) 



{N-2i) 



(13) 



can 



(14) 



Replacing the sum in Eq. [Mlby an integral, using Sterling's approximation for factorials, expand- 
ing the logarithm of the integrand in a Taylor series (up so second order) around the maximum 
which is at Vax = {N/2)\/P'^ — A^/ (1 — P + ^P"^ — A^) and then exponentiating the expansion, 
and finally extending the limits of integration to ±oo (which has a negligible effect on the result 
for > 1) - leads to: 



n(A(t) = 0) 



1 - P + VP2 - A2 



N 



X 



dyexp 



C{1-C)N 



iy ^max) 



where C = VP^ - A^/{1 - P + VP^ - A2). This yields 



2ro 



n (A (t) = 0) 



2ro. 



7rC(l - C)N 



1 - P + VP2 - A2 



-N 



(15) 



(16) 
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where tq is some microscopic time scale. (The factor of 2 in the numerator in Eq. [16] is due to 
the fact that once the motors stop, they have equal probability to move in both directions.) As 
noted before (see section [1]), Eq. [16] predicts an almost exponential dependence of Trev on only 
for constant values of P and A, which was the case in ref. [23]. In the more general case, the 
dependence of t-^^w on N can be derived by calculating the values of P and A as a function of 
and substituting these values into Eq. [THl 

To test the vahdity and accuracy of the analytical expression for Tj-ev, we take the following 
steps: (i) set the model parameters /, a, cui, uj2, tu^, and a to the values used in our computer 
simulations which are given in Table [T] (ii) calculate the probability density patt corresponding to 
these values (Eq. [9]) and use Eqs. [lO] and [H] to calculate P and A over the range of A^ studied in 
the simulations, (iii) substitute the values of P and A into Eq. [161 to obtain function of 

A^, (iv) fit the analytical expression for rrev(A^) to the simulation results plotted in Fig. [2f3. This 
procedure involves two fitting parameters: the microscopic time scale tq appearing in Eq. [16], 
and the group velocity v appearing in the steady-state equation (Eq. [9]). A seemingly reasonable 
choice for the latter would be f = 20 nm/sec, which is where the velocity histogram of the 
bidirectional motion is peaked (see inset of Fig. [3J\). However, the motors slow down before 
each change in their direction of the motion; and because these changes in the directionality are 
fairly rare events, their occurrence probability is likely to be strongly influenced by the short 
periods of slow motion preceding them. Thus, it can be expected that the best fit of Eq. [16] to 
the simulation results is achieved for u < 20 nm/sec. Indeed, for v = 8.2 nm/sec and tq = 680 
msec, we obtain the fitting curve shown in Fig. [6JA., which is an excellent agreement with our 
computational results for the reversal times (plotted in Fig. [2]B and replotted here in Fig. [6]A) 
over the whole range of values of A^ investigated (400 < N < 2400). The steady state probability 
density, Pa.tt{x), on the basis of which t^^v was calculated is shown in Fig. [6l3 for several different 
values of A^ (A^ = 1000 - solid line, A^ = 2000 - dashed line, A^ = 2500 - thick solid line). As can 
be seen from the figure, the detachment rate ui in our simulations is so large that the central 
detachment interval of the unit cell (—a < x < a) is completely depleted of motors. Increasing 
A^ leads to a decrease in the effective attachment rate around the potential maximum, which 
reduces both the number of motors supporting {—1/2 < x < —a) and objecting (a < x < 1/2) 
the motion and leads to the non-monotonic dependence of r on A^. The fitting value of tq = 680 
msec is very close to r* = l/v = 5 nm/(8 nm/sec) = 625 msec, which is the traveling time of the 
motors within a unit cell of the potential (once we set v = 8.2 nm/sec) and, therefore, is also 
the characteristic time scale at which the motors change their "states" (detached, connected and 
supporting the motion, connected and objecting the motion). The remarkable agreement between 
the analytical and simulation results for t^^v should not, however, be allowed to obscure the fact 
that Eq. [TH] is based on a mean field approximation which, in principle, is not suitable for the 
calculating the probabilities of rare fluctuation events (such as velocity reversals in cooperative 
bidirectional movement). The agreement is achieved with effective velocity {v = 8.2 nm/sec) 
which is significantly smaller than the typical velocity measured in the simulations {v = 20 
nm/sec). Therefore, one should not expect the steady state probability density Patt(a^) plotted 
in Fig. [n]B to perfectly match the simulations data. 
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4 Summary 



We use a two-state ratchet model to study the cooperative bidirectional motion of myosin II 
motors on actin tracks with randomly alternating polarities. Our model is an extension of a 
model previously proposed by Badoual et al. to explain the macroscopically large reversal times 
measured in motility assays [23]. These time scales of velocity reversals are orders of magnitude 
longer than the microscopic typical stepping times of individual motors and can be understood as 
a result of collective effects in many-motor systems. The ratchet model that we use assumes that 
the motors are coupled mechanically but act independently, i.e., their binding to and unbinding 
from the cytoskeletal track are statistically uncorrclatcd. These assumptions lead to a predicted 
exponential increase of Trev with A^, the number of motors. Motivated by recent experiments 
which exhibit no such dependence of Trev on A'^ [9], we introduced a modified version of Badoual's 
model which accounts for an additional cooperative effect of the molecular motors and which 
eliminates the exponential increase of Tj-cv with A^. This additional collective effect arises from 
the forces that the motors jointly exert on the actin and the associated elastic energy which 
(within a mean field approximation) scales as E/KbT ~ NNc (where Nc < N is the number 
of attached motors). This scaling relationship implies that the typical energy released when a 
motor is detaching from the track increases linearly with N and, therefore, the detachment rate 
in many-motor systems should be larger than the detachment rate of individual motors. We 
show, both computationally and analytically, that when this effect is taken into account and 
the detachment rate is properly redefined, the characteristic reversal time does not diverge for 
large N. Instead, Trev exhibits a much weaker dependence on N and reaches a maximum at 
intermediate values of N. 

While our model definitely improves the agreement with the experimental results (compared 
to the original model), further improvement is needed in order to eliminate the non-monotonic 
dependence of Trev on N. One step in this direction may be to consider other forms of the off-rate 
cus which are based on more accurate evaluations of the actin elastic energy. In the present work, 
our analysis is based on a mean field approximation which makes the calculation tractable by 
assuming that the detachment rate depends only on (the total number of motors), but not on 
the instantaneous number of attached motors and their locations along the cytoskeletal track. A 
full statistical mechanical treatment is feasible only for small systems, which we plan to report in 
a future publication. As a final remark here we note that the mean field approximation probably 
leads to over-estimation of the effect of the "track-mediated" elastic interactions on the reversal 
times (which may explain the decrease in Trev for large A^). In a non-mean- field calculation the 
motors which release higher energy will detach at higher rates, and the detachment of these 
"energetic" motors will lead to the release of much of the elastic energy stored in the actin track. 
By contrast, in the mean field approximation the contribution of all the connected motors to the 
energy is the same. Therefore, within the mean-field approximation, a larger number of motors 
must be disconnected at a higher frequency, which increases the "stochastic noise" in the system 
that reduced Trev 
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the Joseph and May Winston Foundation Career Development Chair in Chemical Engineering, 
the Israel Cancer Association (grant No. 20070020B) and the Israel Science Foundation (grant 
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Figure 1: point particles (representing the motors) are connected to a rigid rod with equal 
spacing q. The motors interact with the actin track via a periodic, symmetric, saw-tooth potential 
with period / and height H. In each periodic unit, there is a random force of size /ran, pointing 
either to the right or to the left. The motors are subject to these forces only if connected to 
the track. The detachment rate uji is localized in the shaded area of length 2a < I, while the 
attachment rate a'2 is located outside of this region. The off rate U3 is permitted only outside 
the gray shaded area. 




Figure 2: (A) The probability distribution function (PDF) of reversal times corresponding to one 
track realization. The distribution is fitted by a single exponential decay function - see Eq. [H 
(B) The mean reversal time (rrev) as a function of the number of motors N. For each value of N, 
the calculation of (rrev) is based on simulations of 40 different track realizations, where the error 
bars represent the standard deviation of T^ev between realizations. The solid and open circles 
denote the results corresponding to a = 0.002 (our model) and a = (original model presented 
in ref. [23]), respectively. In the latter increases exponentially with (as indicated 

by the solid straight line), while in the former case (xrev) exhibits a non-monotonic behavior (as 
indicated by the dashed line which serves as a guide to the eye) and reaches considerably lower 
values. The triangles denote the results for periodic tracks whose reversal times are always larger 
than those measured on random tracks. 
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Figure 3: The probability distribution functions (PDF) of reversal times corresponding to = 
1000 motors moving on slightly polar tracks. The variable D denotes the difference between the 
fraction of random forces pointing in the favored and disfavored directions. The motion in the 
favored and disfavored directions are analyzed by different PDFs, each of which can be fitted 
by a single exponential form but with distinct reversal times (except for D = where the two 
PDFs coincide with each other). The insets show the velocity histogram corresponding to each 
value of D. As D increases, the histograms become less and less symmetric. 
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Figure 4: The mean reversal times (computed based on simulations of 8 track realizations with 
N = 1000 motors) as a function of D, the difference between the fraction of random forces along 
the track which point in the favored and disfavored directions. The motion in the favored and 
disfavored directions are characterized by the larger (rrev-i) and smaller (rrev-s) reversal times, 
respectively. 
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Figure 5: The steady state probability density, Patt, as a function of the position within a 
unit cell of the periodic potential. The functions plotted in the figure correspond to 2a = 0.76/, 
cug = 0, and (t<;i,u;2) = {v/l,v/l) - thin solid line, {001,002) = {5v/l,5v/l) - dashed line, {001,002) = 
{30v/l,30v/l) - thick solid line. 
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Figure 6: (A) The reversal time (rrev) function of the number of motors A^. The circles denote 
the simulations results (replotted from Fig. [2]). The curve is a fit of the results to Eq. [161 with 
To = 680 msec and v = 8.2 nm/sec. (B) The steady state probability density Pa.tt{x) computed 
for several values of (A^ = 1000 - solid line, A^ = 2000 - dashed line, A^ = 2500 - thick solid 
line). The group velocity of the motors is f = 8.2 nm/sec [as in (A), above], while the model 
parameters are set to the values used in our simulations (see Tabled]). 
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